#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int JacobiRichardson (double ** MA, double * MB, double * MX, double J_ERROR, int J_ORDER, int J_ITE_MAX)
{
    int i, j, iterations = 0;
    double sum, max_error, max_x, aux;
    double * MX2 = (double *) calloc(J_ORDER, sizeof(double));
    double * AUX;
    
    for (i = 0; i < J_ORDER; i++) {
        for (j = 0; j < J_ORDER; j++)
            if (i != j)
                MA[i][j] /= MA[i][i];
        MB[i] /= MA[i][i];
    }
    
    for (i = 0; i < J_ORDER; i++)
        MA[i][i] = 0;
    
    for (i = 0; i < J_ORDER; i++)
        MX[i] = MB[i];
    
    do {
        for (i = 0; i < J_ORDER; i++)
        {
            MX2[i] = MB[i];
            sum = 0.0;
            for (j = 0; j < J_ORDER; j++) sum += MA[i][j] * MX[j];
            MX2[i] -= sum;
        }
        max_error = 0;
        for (i = 0; i < J_ORDER; i++)
        {
            aux = fabs(MX2[i]-MX[i]);
            if (aux > max_error)
                max_error = aux;
        }
        max_x = 0;
        for (i = 0; i < J_ORDER; i++)
        {
            aux = fabs(MX2[i]);
            if (aux > max_x)
                max_x = aux;
        }
        
        AUX = MX2;
        MX2 = MX;
        MX = AUX;
        iterations++;
        
    } while (iterations < J_ITE_MAX && max_error/max_x > J_ERROR);
    
    return iterations;
}

int main (int argc, char ** argv)
{
    /* Declaração de variáveis */
    
    int i, j, iterations;
    int J_ORDER, J_ROW_TEST, J_ITE_MAX;
    double J_ERROR, J_AUX, MBR;
    double ** MA;
    double * MB;
    double * MX;
    double * MR;
    FILE * f;
    
    if (argc != 2)
    {
        printf("Uso: %s nomedoarquivo.txt\n", argv[0]);
        return -1;
    }

    /* Abertura de arquivos */
    f = fopen(argv[1], "r");
    if (!f) return -1;

    /* Leitura dos arquivos */
    fscanf(f, "%d", &J_ORDER);
    fscanf(f, "%d", &J_ROW_TEST);
    fscanf(f, "%lf", &J_ERROR);
    fscanf(f, "%d", &J_ITE_MAX);

    /* Alocação dinâmica e leitura das matrizes */
    MA = (double **) malloc( J_ORDER * sizeof(double *) );
    MX = (double *) malloc( J_ORDER * sizeof(double) );
    MB = (double *) malloc( J_ORDER * sizeof(double) );
    MR = (double *) malloc( J_ORDER * sizeof(double) );

    for (i = 0; i < J_ORDER; i++)
        MA[i] = (double *) malloc(J_ORDER * sizeof(double));

    for (i = 0; i < J_ORDER; i++)    
        for (j = 0; j < J_ORDER; j++)
            fscanf(f, "%lf", &MA[i][j]);
    
    for (i = 0; i < J_ORDER; i++)
        fscanf(f, "%lf", &MB[i]);
    
    for (i = 0; i < J_ORDER; i++)
        MR[i] = MA[J_ROW_TEST][i];
    
    MBR = MB[J_ROW_TEST];
    
    /* Chamada da função Jacobi-Richardson */
    iterations = JacobiRichardson(MA, MB, MX, J_ERROR, J_ORDER, J_ITE_MAX);
    
    /* Teste de linha */
    J_AUX = 0;
    for (i = 0; i < J_ORDER; i++)
        J_AUX += MR[i] * MX[i];

    /* Saída formatada */
    printf ("\n---------------------------------------------------------\n");
    printf ("Iterations: %d\n", iterations);
    printf ("Row test: %d => [%lf] =? %lf", J_ROW_TEST, J_AUX, MBR);
    printf ("\n---------------------------------------------------------\n");

    /* Desalocando memória */
    fclose(f);
    for (i = 0; i < J_ORDER; i++) free(MA[i]);
    free(MA);
    free(MX);
    free(MB);
    free(MR);
    
    return 0;
}
